Two-dimensional Granular Gas of Inelastic Spheres with Multiplicative Driving 
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We study a two-dimensional granular gas of inelastic spheres subject to multiplicative driving pro- 
portional to a power ^(a:)!* of the local particle velocity v(x). The steady state properties of the 
model are examined for different values of <5, and compared with the homogeneous case 5 = 0. 
A driving linearly proportional to v(x) seems to reproduce some experimental observations which 
could not be reproduced by a homogeneous driving. Furthermore, we obtain that the system can be 
homogenized even for strong dissipation, if a driving inversely proportional to the velocity is used. 

PACS: 45.70, 47.50+d, 51.10.+y, 47.11,+j 



Granular materials are a fascinating non-linear, dis- 
sipative, non-equilibrium system Jl|, ||, which are of in- 
terest, first due to their practical importance and second 
because of the challenges they represent theoretically and 
numerically, but experimentally as well. Granular media 
are collections of macroscopic particles with rough sur- 
faces and dissipative, frictional interactions. For the sake 
of simplicity, we restrict ourselves to the smooth situation 
here, and use the inelastic hard sphere (IHS) model |j] [|. 
For more realistic interaction models see Refs. and 
the references therein. 

One of the outstanding effects is the so called clus- 
tering, a self-stabilized density instability due to dissi- 
pation ||,|| where large, dense collections of particles 
coexist with almost empty areas. Clustering occurs in 
initially homogeneous systems and should not 

be confused with the so-called "inelastic collapse" 
the divergence of the collision rate, which is inherent 
to the frequently used hard (rigid) sphere model |3],||, 
Freely cooling systems - mainly examined numeri- 
cally and theoretically (|] - are almost impossible to re- 
alize experimentally. Only recently laboratory experi- 
ments were performed, where clustering could be exam- 
ined in driven systems [|s| — , where also kinetic theory 
approaches complemented by numerical simulations have 
proven to be successful p|,p[[7|,[T^-^5|. 

The driving of a granular material can be realized by 
moving walls || which lead to local heating, or the sys- 
tem can alternatively be driven by a global homogeneous, 
random energy source in different variations [Q, g, [b|, |l5| . 
The latter type of energy input does not exactly resem- 
ble the experimental situation, where a two-dimensional 
(2D) horizontal layer of spheres is agitated by vertical vi- 
brations of the bottom surface and the horizontal degrees 
of freedom are indirectly agitated due to the different 
vertical jump heights of the colliding particles pli|.p^[. 
Thus, the choice of the driving term to put into a the- 
ory for dissipative systems is an open problem and we 
expect that it depends on the nature of the driving (vi- 
brating wall, airflow, Brownian noise, etc.). In order to 



find possible candidates for a realistic energy input, the 
IHS model with an inhomogeneous, multiplicative driving 
is examined in the following. The driving is proportional 
to the local velocity |u(a;)| 5 , with a given power 6 which 
can take both positive and negative values. The classi- 
cal homogeneous driving with S = is contained in our 
approach as a special case. A positive power leads to 
weak energy input for slow particles, e.g. those moving 
collectively inside a cluster ]i"o|,|il|. 

In this letter a system of N three-dimensional spheres 
with radius a and mass m is considered, interacting via 
a hard-core potential and confined to a 2D plane of 
linear extension L, with periodic boundary conditions. 
The degrees of freedom are the positions T"j(t) and the 
translational velocities Vi(t) for each sphere numbered 
by i = 1, . . . , N. The dissipation at a collision is quan- 
tified by a constant normal restitution r. From the mo- 
mentum conservation law we can derive the change of 
linear momentum — (m/2)(l -I- r)vc of particle i which 
collides with particle j, with the normal relative velocity 
Vc 1 ^ — [(vi — Vj) ■ n]n, and the unit vector in normal 
direction n = (r, — Tj)/\ri — rj\. The system is agi- 
tated each time interval At = f7 , with a driving rate 
/dr according to the frequency of the bottom. For ho- 
mogeneous driving, it is costumary to assume /d r S> w, 
where to is the collision frequency of the granular gas 
[^,|^]. This, however, does not correspond to the experi- 
ments, were /d r can be rather small - we will use driving 
frequencies around 100 s" 1 , comparable to experimental 
situations |l(J . Numerical checks with strongly different 
values of /d r lead to a similar behavior of the system even 
for driving frequencies lower than, but of the same order 
as u>, provided that a stationary state is reached. The 
velocity of particle i is modified at each time of agitation 
t so that 



v?(t)=vf(t) + rf\v i (t)\ s v 1 r - s 
v'V(t)=vUt)+rV\v i (t)\ S v r - 5 



(1) 



where the prime on the left hand side indicates the value 
after the driving event. v r is a reference velocity (in 
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this study we use v r = 1ms x ) which allows to de- 
fine the dimensionless translational particle temperature 
T = E/{NT r ), with E = (m/2)X)iIiU- and the refer- 
ence temperature T r = mv 2 . The variance of the uncor- 
rected Gaussian random numbers rf and rf (with zero 
mean) can now be interpreted as a dimensionless driving 
temperature T dr . The stochastic driving rule in Eq. (|l|) 
leads thus to an average rate of change of temperature 



AT/ At = H dl T s with H di . = f dr T dr 



(2) 



The time evolution equation of T was derived for the 
case of a freely cooling granular gas by means of a 
pseudo-Liouville operator formalism [||||. We adopt 
their nomenclature and account for the driving by adding 
Eq. (|f) to the mean field (MF) equation for the transla- 
tional degree of freedom 



—Tit) = -G r AT 3 ' 2 + H dr T s 
at 



(3) 



For our case of a homogeneous monolayer of smooth 

i i 2 

spheres, one has G r — %anyJirT r /mg{v), and A = ~T , 

with the number density n = N/V, the pair corre- 
lation function at contact g{y), and the area fraction 
v = Nna 2 /V covered by particles For the two- 

dimensional layer of spheres we use the approximation 
g{u) = (1 - 7i//16)/(l - vf |J|. For 5 = the driving is 
homogeneous and independent of the local granular tem- 
perature (or velocity). In the case 6 ^ the driving is 
a function of T [0. Imposing 4rT(t) — one gets from 
Eq. (0) the MF temperature in the steady state 



H dl 

gla 



(4) 



the generalization of the Enskog equilibrium solution for 
a driven granular gas |^|, |) . The scaling exponent of T mf 
in Eq. | is 2/3 for 5=0, while it is 2 for 5 = 1. For 5 > 
3/2 the MF theory does not admit a stable equilibrium 
state H[l§. 

The final approach to the steady state can be obtained 
by linearizing Eq. (ph around T mf , what leads to an ex- 
ponentially fast approach to equilibrium with the power 
-[ZAw + 5i/ d r(T mf )' 5_1 ]t, where u = G r ^/T^ '/2 is the 
Enskog collision frequency for elastic particles with tem- 
perature T mf . By inserting Eq. ([!]) in the expression 
for u>, one can express the characteristic relaxation time 
trelax = [• • ■] of T(t) as a function of the model parame- 
ters, which reduces for S = 1 to i TO | ax = (5/2)H dr . Thus, 
for 5 = 1 the characteristic time for the evolution of T to- 
wards its equilibrium value does not depend on A which 
contains all the information about the inelasticity. This 
characteristics is confirmed by numerical simulations. 

Most of our event driven (ED) molecular dynamics 
simulations, see [^|-^) for details, with the driving speci- 
fied in Eq. (|l|) , are first equilibrated without driving and 



with elastic interactions (r = 1), until the velocity distri- 
bution is close to a Maxwellian. Then, dissipation and 
driving are switched on. However, we checked that the 
steady state does not depend on the initial conditions. 
The simulations are performed at fixed volume fraction 
v = 0.34 with N = 1089 (/ dr = 133 s" 1 ) or N = 11025 
C/dr = 67 s" 1 ), and different values of r. In our simula- 
tions we have chosen a = 10~ 3 m and H dl = 1.0 s^ 1 so 
that, for example, G r = %vv r j 1 ^fiiagiv) = 3.1 x 10 3 s~ 1 , 
T mf = 0.0358 and thus w = 2.9 x 10 2 s~\ if r = 0.90 is 
used. With these typical values and a homogeneous driv- 
ing (5 = 0), the model of a driven granular gas is very 
close to a homogeneous state; no clusters are observed 
and the velocity distribution is almost Maxwellian JJ7 1 . 
This is in contradiction to the experimental findings [pLO | , 
and suggests that the correct representation of the driv- 
ing in those experiments is not the homogeneous white 
noise usually implemented. 
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FIG. 1. Snapshots of the particle distribution in the steady 
state for a system of TV = 11025 particles, 8 = 1, v = 0.34, 
and H dl = 1.0s" 1 , with r = 0.97 (a) and r = 0.6 (b). 

In Fig. [I] snapshots of the system's steady state are 
shown for r = 0.97 (a) and r — 0.6 (b), with 5=1. 
Different regimes are observed: A homogeneous state 
exists for very weak dissipation (r = 0.999, not shown 
here), whereas dense, persistent clusters with a crys- 
talline structure, domain boundaries, and vacancies are 
found for higher dissipation (r = 0.97). The region be- 
tween these clusters appears rather homogeneous and di- 
lute (gas-like), very similar to the structures observed in 
experiments [ flO[ . For higher dissipation (r = 0.6) the 
clusters appear less symmetric and are smaller. 
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Note that the clustering in the case S = 1 is qualita- 
tively different from the case of homogeneous driving, and 
it appears already for quite high values of r (see Fig. [l]) . 
Homogeneous driving, in fact, leads to transient clusters, 
i.e. they appear and disappear continuously, while in the 
experiments |lC| and in the case of multiplicative driv- 
ing, the individual clusters are in equilibrium with a gas 
phase and are stable for rather long times. Simulations 
with negative 5, (we used S = —0.5 and 6 = —0.25) give a 
behavior qualitatively similar to the case of homogeneous 
driving (5 = 0). For S = 1, the particles inside the clus- 
ter, with rather small relative velocities, are much less 
agitated than the particles in the surrounding gas phase, 
so that a cluster is stable. For S < 0, the particles in 
the cluster are driven comparatively strong, what leads 
to less stable, dynamic clusters. 

In Fig. H we plot the ratio between the numerical re- 
sults for long times T cq and the theoretical equilibrium 
temperatures T mi as function of r for different S. The 
agreement of the simulations with the MF prediction is 
optimal for r — ► 1. For 5 < 0, the range of agreement 
extends to much smaller r values, i.e. to stronger dissipa- 
tion, as in the case of (5=1 and even in the case 5 = 0. 



function f a (v x ) = fo exp(— B\v x — (v x )\ a ) is used to esti- 
mate the exponent a of the tail of the distribution. The 
results are reported in Fig. [|, where a representative ve- 
locity distribution is shown in the inset. The main out- 
come of our simulations is that the velocity distribution is 
not Gaussian for high inelasticity and that the exponent 
a in the stretched exponential depends on the normal 
restitution r. The non-Gaussian behavior is limited to 
the tail for large r, but seems to extend over the whole 
range of velocities if r is small enough. 
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FIG. 2. Rescaled translational temperature r cq /T mf plot- 
ted against the restitution coefficient r for N = 11025 and 
different values of 5 as given in the insets. Note the different 
axis scaling in (a) and (b). 

We have also performed simulations with 5 = — 1, and 
found that the driving is very singular in the low velocity 
limit, i.e. an excessive amount of energy is given to slow 
particles, creating an inhomogeneous state. 

In the following, we focus on the steady state velocity 
distribution. Therefore, we performed a large series of 
simulations for different values of r, with v = 0.34, 6=1, 
and both N = 1089 and N = 11025. For each value of 
r we performed about 10 3 simulations for N = 1089 and 
about 2 x 10 2 simulations for N = 11025, each with dif- 
ferent initial configurations. The distribution of velocity 
is symmetric and isotropic so that we present only the 
data for the x component. For data analysis, a Gaussian 
with the width obtained from the numerical simulations 
can be superposed to the velocity distribution function 
f(v x ) in order to visualize the differences. Furthermore, 
as a more quantitative approach, a three parameter fit 



FIG. 3. Plot of the exponent a as a function of r, for system 
size N = 1089 (solid triangles) and N = 11025 (open squares). 
In the inset, the dots give the averaged, normalized velocity 
distribution from simulations with N = 1089 and r = 0.9, the 
solid line is a fit according to f a {v x ) with fo = 7.6, B = 18.9 
and a = 1.05(9), and the dashed line is a Gaussian with the 
same standard deviation as the data. 

In Fig. ^|, the exponent a varies continuously between 
0.5 and 2 and, in particular, for N = 1089 and r = 0.9 
we get a = 1.05(9), while for N = 11025 and r = 0.95 
we observe a = 1.02(6). Such an almost exponential 
velocity distribution has recently been observed in ex- 
periments on a vibrated monolayer of spheres, together 
with the cluster structure discussed above |T(J. The size 
dependency of our results, seems reasonable, since also in 
freely cooling systems the value of r below which the in- 
homogeneity and the clustering become important grows 
with the system size Simulations for 5 = —0.25 

and —0.5, show that the velocity distribution is very near 
to a Gaussian even for strong dissipation. 

In order to study the shape of the velocity distribution 
function analytically, we extend the approach used by 
van Noije et al. |5| in the special case 6 = to arbitrary 
6. The Enskog-Boltzmann equation for the freely cool- 
ing gas of spheres, with the standard collision integral 
/(/, /) , has to be extended by the multiplicative driving 
term proportional to the particle velocity. The Enskog- 
Boltzmann equation, see Eq. (26) in 0], corrected by a 
Fokker-Planck diffusion term |l7| , is 



3 



d_ 

di 



f(v 1 ,t)=g(2a)I(f,f) 



2m dv\ 



v?f(vi,t)}. (5) 



Considering the stationary limit df(vi,t)/dt — 0, and 
introducing a scaled distribution function, with the di- 
mensionless velocity c = v/vq and the mean thermal ve- 
locity i>o, one obtains a dimensionless evolution equation 
for the scaled /(«), as in Q. Multiplied by c\ and in- 
tegrated over ci, this leads to a set of equations which 
couple the moment (c p ~ 2+2S ) to the p-th moment fx p of 
the collision term. In the special case 6 = 1 and p = 2, 
one obtains vq = Hdr/igi^^^n). For 5^1 and 5 =/= 
the MF thermal velocity can be obtained by assuming 
(c 25 ) = l(c 2 )] S = < _ 
For high c, one has /(/,/) ~ —P\C\f{ci), with fa = 



rV2 



in 2D and the equation for /(c) becomes 



2d 



d2_ 

dc 2 



1 d 



dc 



c 25 /(c)]=0, (6) 



Inserting a solution of the form f(c) cx exp(— Sc a ) 



we obtain the large c solution 



3-2<5 



and £? 



j^2sy/2dfli/ fi2- ln particular, for 5 = 1 we find a = 1/2 
what corresponds to the exponent we found in the low 
r limit, see Fig. ||. Our explanation is that there is a 
crossover velocity c* above which the tail with exponent 
1/2 appears, which decreases as r decreases. Such a tail 
is not observed numerically for high r values, since c* 
is outside the velocity range we can examine. Simula- 
tions with r = 0.2, N = 1089, and 6 = -0.5, -0.25, and 
0.5 give respectively a = 2.10(9), 1.90(9), and 1.10(7), 
quite near to the values a = 2, 7/4, and 1, given by 
the large c analysis. Since the driving depends on the 
velocity, it acts on all the velocity scales, affecting the 
shape of f(v) in a more complex way than in the case of 
homogeneous driving. Moreover, when inhomogencitics 
are present, the driving will depend on the local ther- 
mal velocity vq(x) and on the local density n(x), while 
homogeneous driving is identical everywhere. Since inho- 
mogeneities are generated by the dissipation, a feed-back 
between driving and dissipation could be at the origin of 
the non- universal behavior of f(v). 

Future research includes a MF and clustering study 
for a wider parameter range, including the case of rough 
particles. Both a more detailed analysis of f(v) and a 
quantitative analysis of the clusters [^0| are in progress. 
A microscopic justification for the multiplicative driving 
is still lacking and requires more detailed experimental 
or three-dimensional simulation |lj] studies. However, 
we suggest a possible experimental check of our ideas. 
From MF theory the scaling of the equilibrium tempera- 
ture of the vertically vibrated gas with the area fraction 



v is given by T n 



Experimental mea- 



surement of the equilibrium temperature of the vertically 
vibrated monolayer, for different values of v, could allow 
to estimate the value of 6 and to verify the hypothesis of 



a multiplicative effective driving. 
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